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Abstract. We present a probabilistic generative approach for construct- 
ing topographic maps of light curves from eclipsing binary stars. The 
model defines a low-dimensional manifold of local noise models induced 
by a smooth non-linear mapping from a low-dimensional latent space into 
the space of probabilistic models of the observed light curves. The local 
noise models are physical models that describe how such light curves 
are generated. Due to the principled probabilistic nature of the model, 
a cost function arises naturally and the model parameters are fitted via 
MAP estimation using the Expectation-Maximisation algorithm. Once 
the model has been trained, each light curve may be projected to the 
latent space as the the mean posterior probability over the local noise 
models. We demonstrate our approach on a dataset of artificially gener- 
ated light curves and on a dataset comprised of light curves from real 
observations. 
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1 Introduction 

The Generative Topographic Map algorithm (GTM) has been introduced as 
a probabilistic analog to SOM [2], seeking to address certain of its Umitations 
such as the absence of a cost function. The GTM formulates a mixture of spher- 
ical Gaussians densities constrained on a smooth image of a low-dimensional 
latent space. Each point in the latent space is mapped via a smooth non-linear 
mapping to its image in the high-dimensional data space. This image plays the 
role of the mean of a local spherical Gaussian noise model that is responsible 
for modelling the density of data points in its vicinity. The GTM can be readily 
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extended to structured data by adopting alternative formulations of noise mod- 
els in the place of Gaussian densities. Such extensions have been proposed in 
[3] for the visualisation of symbolic sequences and in 4 for the visualisation of 
tree-structured data. 

Here we present a further extension of the GTM to a novel data type, namely 
light curves that originate from eclipsing binary systems. Binary stars are gravi- 
tationally bound pairs of stars that orbit a common centre of mass. Astronomical 
observations suggest that almost half of the stars are binary ones. Thus, study- 
ing such systems procures knowledge for a significant proportion of stars. Binary 
stars are important to astrophysics because they allow calculation of fundamen- 
tal quantities such as masses and radii, and are important for the verification 
of theoretical models for stellar formation and evolution. A particular subclass 
of binary stars are eclipsing binary stars. The luminosity of such stars varies 
over time and forms a graph called light curve. Light curves are important be- 
cause they provide information on the characteristics of stars and help in the 
identification of their type. 

2 Physical Model For Eclipsing Binaries 

The physical model that generates light curves from eclipsing binary systems 
is described by the following set of parameters: mass Mi e [0.5, 100] (in solar 
masses) of the primary star (star with highest mass of the pair), mass ratio 
q G [0,1] (hence mass of secondary star is AI2 = qMi), eccentricity e G [0,1] 
of the orbit and period p e [0.5, 100] measured in days, all of which specify the 
shape of the orbit. Furthermore, two angles describing the orientation of the 
system are necessary [5] which are known as the inclination i 6 [0, ^] and the 
argument of periastron uj S [0, 2tt] (see Fig. [T|). Inclination describes the angle 
between the plane of the sky and the orbital plane and periastron is the angle 
Lo G [0, 27r] that orients the major axis of the elliptic orbit within its plane, 
that is Lu is measured within the orbital plane. Finally, a third angle known 
as the longitude of ascending node (]7 £ [0, Stt]) is necessary for the complete 
description of a binary system. However, since it has no effect on the observed 
light curves, we omit it from the model. We collectively denote these parameters 
by vector 0. 

The mass M of each star relates to the luminosity L radiated by a surface 
element [6] of the star according to L = M^-^ . Moreover, masses relate to the 
radii R of the stars via: 



These relations show that the primary star is the most luminous one and the one 
with the greatest area of the pair (a star appears as a disc to an observer). Thus, 
the observed area of a star is A = -kB? and the observed luminosity is LttE? . 
Henceforth, we index quantities related to the primary star by 1 (e.g. primary 
mass is Mi) and 2 for the secondary star. 



R = 



{ 



200.053+0.977 logio(Af)^ ^ 1.728; 

200.153+0.556 logio(M)^ Otherwise. 



(1) 
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It is shown from Newton's laws that the orbits of an object in the gravita- 
tional field of another object is a conic section of eccentricity e. Here we are 
interested in the case where < e < 1 that corresponds to closed orbits. We 
formulate two-body systems as systems where one body is fixed and the other 
is in orbital motiory. 

The position of the orbiting body is calculated by Kepler's equation as the 
distance r from the fixed companion star on the elliptical orbit [5], 

where t is time and a is the semi-major axis of the ellipse calculated by Kepler's 
third law. Point 77 in Fig. [T] is the periastron, the point where the distance 
between the orbiting and fixed body is minimal. Angle 9 is the angle between 
the radius and the periastron. Knowledge of 9 would allow us to determine the 
position of the orbiting body. Angle 9 is indirectly inferred via an auxiliary circle 
centered at the center of the ellipse O and radius equal to semi-major axis. Point 
Q is the vertical projection of the orbiting body's position P to the auxiliary 
circle. Angle E is called the eccentric anomaly and is given by Kepler's equation 

U- 

E{t) = esmE{t) + —{t- t), (3) 

where r is the instance of time that the body was at the periastron. Kepler's 
equation does not admit an analytical solution but can be approximated through 
the Newton-Raphson method. By geometrical arguments it is shown that the 
relation between the true and eccentric anomaly reads: 

tan^ = [(l + e)/(l-e)]Uan(^) (4) 

By knowledge of 9 we can determine the position of the second star on the orbit 
using ^ and Q. These positions correspond to the orbital plane and must be 
projected to the plane of the observer in the form of Cartesian coordinates [5]: 

X{t) ^ r(t)(cos(J7) cos(w -I- 9{t)) - sm{f2) sin(w + 9{t)) cos(i)), (5) 
Y{t) = r{t){cos{[2) cos(w -I- 9{t)) + cos{f2) sin(w -I- 9{t)) cos(«)), (6) 
Z(t) = r(t) sin(w -|- 9{t)) sin(z), (7) 

which concludes the determination H of positions of the stars with respect to the 
observer. 



* It is shown in [B] that in the relative motion system, the eccentricity, period and semi- 
major of the moving body's orbit are equal to their counterparts in the two-body 
system, and only the masses transform. 

^ The angle J7 does influence the position of the orbiting body. However, it does not 
have an influence on the light curve and thus we treat it as a constant = 0. 
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Fig. 1. Angles orientating the orbital plane with respect to the plane of sky, 
and angles associated with the orbits. Adapted from [5]. 



An observer of the binary system receives a variable luminosity from the 
eclipsing binary system that plotted against time forms a light curve. This vari- 
ability is due to the eclipses that occur when one body passes in front (in the 
line of sight of the observer) of the other. This is illustrated in Fig. [D When no 
eclipse occurs (positions a,g) the luminosity is equal to the sum of the luminosi- 
ties radiated from the two bodies. The curved parts of the light curve occur due 
to partial occlusions. Two eclipses take place at each period, one primary eclipse 
(position d), when the most luminous body of the pair is obscured the most, 
and a secondary eclipse (position j), when the most luminous body obscures its 
companion the most. 

Obscured parts of the disks of the stars can be calculated via geometrical 
arguments. The obscured area of each star at time t is denoted by AAi{t) and 
AA2 (t) . The luminosity fg {t) received by the observer at time t depends on the 
luminosities Li, areas Ai and obscured area^ AAi via 

fg{t) = LMi - ^Ai{t)) + L2{A2 - AA2{t)). (8) 



® see http://www.physics.sfasu.edu/astro/ebstar/ebstar.html Last access on 12-0-07. 
Recall that i — 1 and i = 2 index the primary and secondary stars, respectively 
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(R1 + R2|> > |R1^2 -R2"2r[l/2] (Rl + R2) >A >(RI'>2 - R2'>2)^(l/2) 



Fig. 2. Positions of stars (relative to observer's line of sight) and corresponding 
light curve phases. 

3 Noise Model for Light Curves 

Based on the physical model a probabilistic generative noise model arises natu- 
rally. Observed light curves, denoted by O, are noisy signals: 

0{t)^ fgit)+e{t), (9) 

where e is i.i.d. Gaussian noise with variance a^. Thus, we regard a light curve 
O of period p(0) sampled at times t £ T = {ti = 0,t2,---,iT = p(0)} as a 
realisation drawn from a multivariate spherical normal distribution. We denote 
the noise model associated with parameters by p{0\f{.; 6), cr^) or simply by 

p{o\e). 

4 Model for Topographic Organisation 

The starting point of our model formulation is the form of a mixture model 
composed of C noise models as described in section [3l 

c 

p(O|0) =^P(c)p(O|0,), (10) 

c=l 

where P(c) are the mixing coefficients, @ encapsulates all parameter vectors 
{dc}c=i:C and p{O\0c) corresponds to the c— th model component with param- 
eter vector Be- We simplify notation p{0\6c) to p(0|c). Assuming that dataset 
V contains N independently generated fluxes O^"-*, the posterior of the is 
expressed as: 

N N C 

p{0\v) a p{&) n \0) = p{0) n E ic) (11) 

n— 1 n— 1 c— 1 
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where the mixing coefficients can be ignored as P{c) = ^. 

Topographic organisation is introduced in the spirit of the GTM [J] by re- 
quiring that the component parameter vectors 6c correspond to a regular grid 
of points Xc,c = 1,...,C, in the two dimensional latent space V = [— f , 1]^. 
A smooth nonlinear function F maps each point a; G V to a point r{x) that 
addresses a model p{-\x). Points r{x) are constrained on a two-dimensional 
manifold A4 that is embedded in space 7i, the space of parametrisations of our 
noise models. Since the neighbourhood of .T-images of x is preserved due to con- 
tinuity of r, a topographic organisation emerges for the models p{-\x). Function 
7^ is realised as a RBF network [T]: 

r{x) = W4>{x), (12) 

where matrix W E R^^^ contains the free parameters of the model (6 is the 
number of parameters in {Mi,q,e,i,uj, p}), and </>(.) = {(t)i{.), ...,(f>K{-))'^ ,<pk{-) '■ 
^ R is an ordered set of K nonlinear smooth basis functions. However, this 
mapping may produce invalid parameter vectors, since the output of the RBF 
network is unbounded. We therefore redefine mapping 7^ as: 

r{x)=Ag{W<l>{x))+v, (13) 

where: 



— g a vector- valued version of the sigmoid function that "squashes" each ele- 
ment in [0, 1]: 



[l + exp(-?/i) ' l-|-exp(-y2)' " ' ' l-f exp(-yY) 

— j4 is a diagonal matrix that scales parameters to the appropriate range. A has 
as diagonal elements the length of range (0™"^ — gi™™) fQj- each parameter, 
sothat^ = dmg((100-0.5), (1-0), (1-0), (27r-0), (f-0), (100-0.5)). 

— vector V shifts the parameters to the appropriate interval, v contains the 
minimum value 6*™*" for each parameter 0^: v = [0.5, 0, 0, 0, 0, 0.5]'^. 

The redefined mapping F now takes a point x in space V to a valid parameter 
vector r{x) that addresses a noise model in Ai. Thus, has become a function 
of the weight matrix W of the RBF network, 0{W). Hence, the logarithm of 
the posterior from (jlip now reads: 

N C 

\ogp{e{W)\V) cxlogp(0(W^)) + ^log^p(O(")|a;c). (15) 

n— 1 c— 1 

Figure [3] summarises the model formulation. Each point x of the visualisation 
space V is non-linearly and smoothly mapped via F to model parameters that 
identify the corresponding noise model p{-\x). These parameters are constrained 
on a two-dimensional manifold embedded in the space of all possible 



(14) 
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parametrisations of our noise model. In the spirit of [T], the model can be used 
to visualise observed fluxes O by calculating the posterior probability of each 
grid point Xc G V, given O: 



p{Xc\0) 



P{x,)p{0\x,) 
p{0) 



Pix,)p{0\x,) 



p{0\Xc) 



Ec'=iPMp{0\x,,) j:%iPiO\xc')' 



(16) 



Each observed flux O is then represented in the visualisation space V by a 
point proj{0) G V given by the expectation of the posterior distribution over 
the grid points: 



proj{0) = ^p{xc\0)xc. 



(17) 










v 





~1 ' coordinate vector 

Latent space [Xp 
(computer screen) 



Distribution space 




Fig. 3. Formulation of the topographic mapping model. 



We train our model in the MAP estimation framework with a physically mo- 
tivated prior p{0) obtained from relevant literature |7|8|9|10j . To that purpose 
we employ the EM algorithm. Note that, due to the nature of the physical model 
formulation in sections [2] and [31 the M-step cannot be carried out analytically, 
nor can the derivatives of expected complete-data log-posterior with respect to 
the RBF network parameters W be analytically obtained. However, the EM al- 
gorithm does not necessarily require that an optimum is achieved in the M-step; 
it is sufficient that the likelihood is merely improved |llj . For our purposes we 
resort to numerical optimisation by employing a (H-l) evolutionary strategy de- 
scribed in [12 . The fitness function for the evolutionary strategy is the expected 
complete-data log-posterior. 
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5 Experiments 
5.1 Datasets 

We performed experiments on two datasets. Dataset J is a synthetic dataset 
that consists of 200 hght curves (fluxes). A common set of model parameters, 
{Ml — 5,q = 0.8, e = 0.3, i = was defined. However, two distinct values 
pi = 2,p2 = 5 of period and oji — Q,uj2 — |7r of argument of periastron were 
used, to create 4 classes of light curves (50 in each class) by the combinations 
of these values, {pi,p2} x {uji,uj2}- The discerning characteristic of each class 
is the position of each secondary eclipse and the widths of the eclipses. Each 
light curve was then generated from these four "prototypical" parameter settings 
corrupted by a Gaussian noise. Gaussian noise was also subsequently added to 
the generated light curves to simulate observational errors. 

Dataset 2 consists of light curves from real observations obtained from two 
resources availabl^on the WWW: the Catalogue and Archive of Eclipsing Bina- 
ries at http://ebola.eastern.edu/land the All Sky Automated Survey. Dataset 2 



was preprocessed before training using local linear interpolations. Preprocessing 
is necessary as one needs to account for gaps in the monitoring process and for 
overlapping observations. Light curves must also be phase-shifted so that their 
first point is the primary eclipse and resampled to equal length as described in 
section 131 Finally, the light curves were resampled at T = 100 regular intervals 
which was judged an adequate sample rate. 



5.2 Training 

The lattice was a 10 x 10 regular grid (i.e. C = 100) and the RBF network 
consisted of M = 17 basis functions; 16 of them were Gaussian radial basis 
functions of variance = 1 centred on a 4 x 4 regular grid in V = [0, 1]^, and 
one was a bias term. The variance of the observation noise in the local models 
p{0\x) was set to ct^ = 0.075. 



5.3 Results 

Fig. m presents the topographic map constructed for the synthetic dataset. Each 
point stands for a light curve projected to latent visualisation space V and is 
coloured according to class membership. The class memberships of synthetic 
fluxes were not used during the training process. Also, next to each cluster, 
a typical light curve has been plotted. The classes have been identified and 
organised appropriately, each occupying one of the four corners of the plot. 

Fig. [5] presents the topographic map constructed for the dataset of real ob- 
served light curves. The red curves are the data projected against the underlying 
local noise models displayed in black. Several interesting observations can be 
made about the topographic formation of the light curves on the resulting map. 



Last accessed on the 12th September 2007. 
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Fig. 4. Visualisation of synthetic dataset. A representative light curve is plotted 
next to each cluster. 

In the lower right-hand corner binary systems of large periods are found. The 
median period of the systems in our sample is 2.7 days, and binaries like V459 
Gas, with a period of 8.45 days lie in this corner. Systems with short period 
have the appearance of a wide V-shaped eclipse in the shape of their light curve, 
and inhabit the top and left edges of the map, e.g. WY Hya (Period: 0.7 days) 
and RT And (Period: 0.6 days). At the lower left of the map, we find systems 
with high eccentricity, e.g. V1647 Sgr. High eccentricity causes the light curve 
to appear assymetric, so that the period of the eclipse occurs further and fur- 
ther away from the center. On the other hand, very symmetric curves indicate 
orbits of low eccentricity (more circular) and low mass-ratio (stars of similar 
mass), and indeed we find systems like DM Vir (e — 0.03, mass ratio=l) and 
CD Tau (e — 0.0, mass ratio=:l.G5) in the cluster in the lower-right hand corner 
of the map. Finally, low-inclination systems, occupy the top left-hand corner of 
the map, and these orbits will have very shallow eclipses as the companion star 
barely eclipses the primary star. 

6 Conclusions 

We have presented a model-based probabilistic approach for the visualisation of 
eclipsing binary systems. The model is formulated as a constrained-mixture of 
physically motivated noise models. As a consequence, a clear cost function nat- 
urally arises which drives the optimisation of the model. In our experiments we 
have demonstrated that the resulting maps can be interpreted in a transparent 
way by inspecting the underlying local noise models. Furthermore, modification 
and refinement of the local noise models is possible, to account for greater phys- 
ical fidelity by incorporating physical aspects for non-spherical stars and even 
more sophisticated phenomena such as gravity darkening. 
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Fig. 5. Visualisation of dataset 2 of real data. Light curves in red are the 
projected real data and light curves in black are the light curves of the underlying 
local noise models. 
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